module ffsl_mod

	use const_mod
	use namelist_mod
	use mesh_mod
	use state_mod
	use tracer_mod
	use parallel_mod
	use operators_mod

	implicit none

	private
	public ffsl_init
	public calc_tracer_mass_flux_hori_ffsl
	public calc_tracer_mass_flux_vert_ffsl

	interface
		subroutine calc_hori_flux_interface(tracer, mx, my, mfx, mfy, dt)
			import tracer_type, r8, mesh
			type(tracer_type), intent(in):: tracer
			real(r8), intent(in) ::   mx(mesh%full_lon_lb:mesh%full_lon_ub, &
				 	 			  								 mesh%full_lat_lb:mesh%full_lat_ub, &
				 	 			  								 mesh%full_lev_lb:mesh%full_lev_ub)
	    real(r8), intent(in) ::   my(mesh%full_lon_lb:mesh%full_lon_ub, &
				 	 			  								 mesh%full_lat_lb:mesh%full_lat_ub, &
				 	 			  								 mesh%full_lev_lb:mesh%full_lev_ub)
			real(r8), intent(out) :: mfx(mesh%half_lon_lb:mesh%half_lon_ub, &
																	 mesh%full_lat_lb:mesh%full_lat_ub, &
																	 mesh%full_lev_lb:mesh%full_lev_ub)
			real(r8), intent(out) :: mfy(mesh%full_lon_lb:mesh%full_lon_ub, &
																	 mesh%half_lat_lb:mesh%half_lat_ub, &
																	 mesh%full_lev_lb:mesh%full_lev_ub)
			real(8), intent(in) :: dt
		end subroutine calc_hori_flux_interface

		pure real(r8) function slope_interface(f)
			import r8
			real(r8), intent(in) :: f(-1:1)
		end function slope_interface
	end interface

	procedure(calc_hori_flux_interface), pointer :: calc_hori_flux => null()
	procedure(slope_interface         ), pointer :: slope          => null()

contains

	subroutine ffsl_init()

		select case(ffsl_flux_type)
		case ('van_leer')
			calc_hori_flux => calc_hori_flux_van_leer
		case ('ppm')
			calc_hori_flux => calc_hori_flux_ppm
		case default
			print*, "Unknown flux_type " // trim(ffsl_flux_type) // '!'
			stop
		end select

		select case (ffsl_limiter_type)
		case ('mono')
			slope => slope_mono
		case ('pd')
			slope => slope_pd
		case ('simple')
			slope => slope_simple
		case default
			print*, "Unknown limiter_type " // trim(ffsl_limiter_type) // '!'
			stop
		end select

	end subroutine ffsl_init

	subroutine calc_tracer_mass_flux_hori_ffsl(old_tracer, old_state, old_q,  q, qmfx, qmfy, dt, do_limiter)

		type(tracer_type), intent(inout) :: old_tracer
		type(state_type), intent(in) :: old_state
		real(r8), intent(in) :: old_q(mesh%full_lon_lb:mesh%full_lon_ub, & ! unused dummy argument
														      mesh%full_lat_lb:mesh%full_lat_ub, &
														      mesh%full_lev_lb:mesh%full_lev_ub)
		real(r8), intent(in) :: q(mesh%full_lon_lb:mesh%full_lon_ub, &
														  mesh%full_lat_lb:mesh%full_lat_ub, &
														  mesh%full_lev_lb:mesh%full_lev_ub)
		real(r8), intent(out) :: qmfx(mesh%half_lon_lb:mesh%half_lon_ub, &
																  mesh%full_lat_lb:mesh%full_lat_ub, &
																  mesh%full_lev_lb:mesh%full_lev_ub)
		real(r8), intent(out) :: qmfy(mesh%full_lon_lb:mesh%full_lon_ub, &
																  mesh%half_lat_lb:mesh%half_lat_ub, &
																  mesh%full_lev_lb:mesh%full_lev_ub)
		real(8), intent(in) :: dt
		logical, intent(in), optional :: do_limiter ! unused dummy argument
		integer i, j, k 
    real(r8) pole

		call calc_hori_cfl(old_tracer, dt)
    call calc_hori_div(old_tracer, dt)

		associate (divx => old_tracer%divx, &
						   divy => old_tracer%divy, &
						   qx   => old_tracer%qx  , &
						   qy   => old_tracer%qy)
		call calc_hori_flux(old_tracer, q, q, qmfx, qmfy, dt)
		do k = mesh%full_lev_ibeg, mesh%full_lev_iend
			do j = mesh%full_lat_ibeg + 1, mesh%full_lat_iend - 1
				do i = mesh%full_lon_ibeg, mesh%full_lon_iend
					qx(i,j,k) = q(i,j,k) - 0.5 * ((qmfx(i,j,k) - qmfx(i-1,j,k)) / mesh%de_lon(j) - q(i,j,k) * divx(i,j,k)) * dt
					qy(i,j,k) = q(i,j,k) - 0.5 * ((qmfy(i,j  ,k) * mesh%le_lat(j  ) - &
						                             qmfy(i,j-1,k) * mesh%le_lat(j-1)) / mesh%area_cell(j) - q(i,j,k) * divy(i,j,k)) * dt
				end do
			end do
		end do
		! south pole
		do k = mesh%full_lev_ibeg, mesh%full_lev_iend
			j = mesh%full_lat_ibeg
			pole = sum(qmfy(mesh%full_lon_ibeg:mesh%full_lon_iend,j,k)) * mesh%le_lat(j) / mesh%num_full_lon / mesh%area_cell(j)
			do i = mesh%full_lon_ibeg, mesh%full_lon_iend
				qx(i,j,k) = q(i,j,k)
				qy(i,j,k) = q(i,j,k) - 0.5_r8 * (pole - divy(i,j,k) * q(i,j,k)) * dt
			end do
		end do
		! north pole
		do k = mesh%full_lev_ibeg, mesh%full_lev_iend
			j = mesh%full_lat_iend
			pole = sum(-qmfy(mesh%full_lon_ibeg:mesh%full_lon_iend,j-1,k)) * mesh%le_lat(j-1) / mesh%num_full_lon / mesh%area_cell(j)
			do i = mesh%full_lon_ibeg, mesh%full_lon_iend
				qx(i,j,k) = q(i,j,k)
				qy(i,j,k) = q(i,j,k) - 0.5_r8 * (pole - divy(i,j,k) * q(i,j,k)) * dt
			end do
		end do
		call fill_merid_halo_cell(qx, all_halo=.true.)
		call fill_zonal_halo_cell(qy, all_halo=.true.)
		! Run outer flux operator
		call calc_hori_flux(old_tracer, qy, qx, qmfx, qmfy, dt)
		end associate

	end subroutine calc_tracer_mass_flux_hori_ffsl

	subroutine calc_tracer_mass_flux_vert_ffsl(old_tracer, old_state, old_q, q, we, mfz, dt, do_limiter)

		type(tracer_type), intent(inout) :: old_tracer
		type(state_type), intent(inout) :: old_state
		real(r8), intent(in) :: old_q(mesh%full_lon_lb:mesh%full_lon_ub, &
			 	 			  								  mesh%full_lat_lb:mesh%full_lat_ub, &
			 	 			  								  mesh%full_lev_lb:mesh%full_lev_ub)
		real(r8), intent(inout) :: q(mesh%full_lon_lb:mesh%full_lon_ub, &
			 	 			  								 mesh%full_lat_lb:mesh%full_lat_ub, &
			 	 			  								 mesh%full_lev_lb:mesh%full_lev_ub)
		real(r8), intent(in) :: we(mesh%full_lon_lb:mesh%full_lon_ub, &
															 mesh%full_lat_lb:mesh%full_lat_ub, &
															 mesh%half_lev_lb:mesh%half_lev_ub)
		real(r8), intent(out) :: mfz(mesh%full_lon_lb:mesh%full_lon_ub, &
																 mesh%full_lat_lb:mesh%full_lat_ub, &
																 mesh%half_lev_lb:mesh%half_lev_ub)
		real(8), intent(in) :: dt
		logical, intent(in), optional :: do_limiter ! unused dummy argument
		integer i, j, k, ku, ci
		real(r8) cf, dq, s1, s2, s3, ds1, ds2, ds3, ql, q6
    
    call calc_vert_cfl(we, old_state%m_lev, old_tracer%cflz, dt)
    associate (cflz => old_tracer%cflz)
		select case (ffsl_flux_type)
		case ('van_leer')
			k = mesh%half_lev_ibeg + 1
			do j = mesh%full_lat_ibeg, mesh%full_lat_iend
				do i = mesh%full_lon_ibeg, mesh%full_lon_iend
					if (cflz(i,j,k) < 0) then! van_leer
						ku = k 
						dq = slope(q(i,j,ku-1:ku+1))
						mfz(i,j,k) = cflz(i,j,k) * (q(i,j,ku) - dq * 0.5_r8 * (1 + cflz(i,j,k))) * old_state%m_lev(i,j,k)
					else if (cflz(i,j,k) > 0) then ! upwind1
						ku = k - 1
						mfz(i,j,k) = cflz(i,j,k) * q(i,j,ku) * old_state%m_lev(i,j,k)
					else 
						mfz(i,j,k) = 0
					end if
				end do
			end do

			k = mesh%half_lev_iend - 1
			do j = mesh%full_lat_ibeg, mesh%full_lat_iend
				do i = mesh%full_lon_ibeg, mesh%full_lon_iend
					if (cflz(i,j,k) < 0) then ! upwind1
						ku = k
						mfz(i,j,k) = cflz(i,j,k) * q(i,j,ku) * old_state%m_lev(i,j,k)
					else if (cflz(i,j,k) > 0) then ! van_leer
						ku = k - 1
						dq = slope(q(i,j,ku-1:ku+1))
						mfz(i,j,k) = cflz(i,j,k) * (q(i,j,ku) + dq * 0.5_r8 * (1 - cflz(i,j,k))) * old_state%m_lev(i,j,k)
					else 
						mfz(i,j,k) = 0
					end if
				end do
			end do

			do k = mesh%half_lev_ibeg + 2, mesh%half_lev_iend - 2
				do j = mesh%full_lat_ibeg, mesh%full_lat_iend
					do i = mesh%full_lon_ibeg, mesh%full_lon_iend
						ci = int(cflz(i,j,k))
						cf = cflz(i,j,k) - ci
						if (cflz(i,j,k) > 0) then
							ku = k - ci - 1
							dq= slope(q(i,j,ku-1:ku+1))
							mfz(i,j,k) = cf * (q(i,j,ku) + dq * 0.5_r8 * (1 - cf)) * old_state%m_lev(i,j,k)
						else if (cflz(i,j,k) < 0) then
							ku = k - ci
							dq = slope(q(i,j,ku-1:ku+1))
							mfz(i,j,k) = cf * (q(i,j,ku) - dq * 0.5_r8 * (1 + cf)) * old_state%m_lev(i,j,k)
						else 
							mfz(i,j,k) = 0
						end if
					end do
				end do
			end do
		case ('ppm')
			! Set upper and lower boundary conditoins: solid wall
			do k = 1, mesh%lev_halo_width
				q(:,:,mesh%full_lev_ibeg - k) = q(:,:,mesh%full_lev_ibeg) ! + 1 - k)
			end do
			do k = 1, mesh%lev_halo_width
				q(:,:,mesh%full_lev_iend + k) = q(:,:,mesh%full_lev_iend) ! + 1 - k)
			end do
			do k = mesh%half_lev_ibeg + 1, mesh%half_lev_iend - 1
				do j = mesh%full_lat_ibeg, mesh%full_lat_iend
					do i = mesh%full_lon_ibeg, mesh%full_lon_iend
						if (cflz(i,j,k) > 0) then
							ku = k - 1
							s1 = 1.0 - cflz(i,j,k)
							s2 = 1.0
							ds1 = s2    - s1
						  ds2 = s2**2 - s1**2
					    ds3 = s2**3 - s1**3
					    call ppm(q(i,j,ku-2:ku+2), ql, dq, q6)
					    mfz(i,j,k) = (ql * ds1 + 0.5 * dq * ds2 + q6 * (0.5 * ds2 - ds3 / 3.0)) * old_state%m_lev(i,j,k)
					  else if (cflz(i,j,k) < 0) then
							ku = k 
							s1 = 0
							s2 = -cflz(i,j,k)
						  ds1 = s2    - s1
					    ds2 = s2**2 - s1**2
					    ds3 = s2**3 - s1**3
					    call ppm(q(i,j,ku-2:ku+2), ql, dq, q6)
					    mfz(i,j,k) = -(ql * ds1 + 0.5 * dq * ds2 + q6 * (0.5 * ds2 - ds3 / 3.0)) * old_state%m_lev(i,j,k)
						else
							mfz(i,j,k) = 0
						end if
					end do
				end do
			end do
		! case ('mix')
		!   k = mesh%half_lev_ibeg + 1
		! 	do j = mesh%full_lat_ibeg, mesh%full_lat_iend
		! 		do i = mesh%full_lon_ibeg, mesh%full_lon_iend
		! 			if (cflz(i,j,k) > 0) then ! upwind1
		! 				ku = k - 1
		! 			  mfz(i,j,k) = cflz(i,j,k) * q(i,j,ku)
		! 			else if (cflz(i,j,k) < 0) then ! van_leer
		! 				ku = k
		! 			  dq = slope(q(i,j,ku-1:ku+1))
		! 			  mfz(i,j,k) = cflz(i,j,k) * (q(i,j,ku) - dq * 0.5_r8 * (1 + cflz(i,j,k)))
		! 			else
		! 				mfz(i,j,k) = 0
		! 			end if
		! 		end do
		! 	end do

		!   k = mesh%half_lev_ibeg + 2
		! 	do j = mesh%full_lat_ibeg, mesh%full_lat_iend
		! 		do i = mesh%full_lon_ibeg, mesh%full_lon_iend
		! 			if (cflz(i,j,k) > 0) then ! van_leer
		! 				ku = k - 1
		! 				dq = slope(q(i,j,ku-1:ku+1))
		! 				mfz(i,j,k) = cflz(i,j,k) * (q(i,j,ku) + dq * 0.5_r8 * (1 - cflz(i,j,k)))
		! 			else if (cflz(i,j,k) < 0) then ! ppm
		! 				ku = k
		! 				s1 = 0
		! 				s2 = -cflz(i,j,k)
		! 				ds1 = s2    - s1
		! 				ds2 = s2**2 - s1**2
		! 				ds3 = s2**3 - s1**3
		! 				call ppm(q(i,j,ku-2:ku+2), ql, dq, q6)
		! 				mfz(i,j,k) = -(ql * ds1 + 0.5 * dq * ds2 + q6 * (0.5 * ds2 - ds3 / 3.0))
		! 			else 
		! 				mfz(i,j,k) = 0
		! 			end if
		! 		end do
		! 	end do

		!   k = mesh%half_lev_iend - 2
		! 	do j = mesh%full_lat_ibeg, mesh%full_lat_iend
		! 		do i = mesh%full_lon_ibeg, mesh%full_lon_iend
		! 			if (cflz(i,j,k) > 0) then ! ppm 
		! 				ku = k - 1
		! 				s1 = 1.0 - cflz(i,j,k)
		! 			  s2 = 1.0
		! 			  ds1 = s2    - s1
		! 		    ds2 = s2**2 - s1**2
		! 	      ds3 = s2**3 - s1**3
		! 	      call ppm(q(i,j,ku-2:ku+2), ql, dq, q6)
		! 	      mfz(i,j,k) = ql * ds1 + 0.5 * dq * ds2 + q6 * (0.5 * ds2 - ds3 / 3.0)
		! 			else if (cflz(i,j,k) < 0) then ! van_leer
		! 				ku = k
		! 				dq = slope(q(i,j,ku-1:ku+1))
		! 				mfz(i,j,k) = cflz(i,j,k) * (q(i,j,ku) - dq * 0.5_r8 * (1 + cflz(i,j,k)))
		! 			else
		! 				mfz(i,j,k) = 0
		! 			end if
		! 		end do
		! 	end do
				
		!   k = mesh%half_lev_iend - 1
		! 	do j = mesh%full_lat_ibeg, mesh%full_lat_iend
		! 		do i = mesh%full_lon_ibeg, mesh%full_lon_iend
		! 			if (cflz(i,j,k) > 0) then ! van_leer
		! 				ku = k - 1
		! 				dq = slope(q(i,j,ku-1:ku+1))
		! 				mfz(i,j,k) = cflz(i,j,k) * (q(i,j,ku) + dq * 0.5_r8 * (1 - cflz(i,j,k)))
		! 			else if (cflz(i,j,k) < 0) then ! upwind1
		! 				ku = k
		! 				mfz(i,j,k) = cflz(i,j,k) * q(i,j,ku)
		! 			else
		! 				mfz(i,j,k) = 0
		! 			end if
		! 		end do
		! 	end do
		! 		 ! ppm
		! 	do k = mesh%half_lev_ibeg + 3, mesh%half_lev_iend - 3
		! 		do j = mesh%full_lat_ibeg, mesh%full_lat_iend
		! 			do i = mesh%full_lon_ibeg, mesh%full_lon_iend
		! 				if (cflz(i,j,k) > 0) then
		! 					ku = k - 1
		! 					s1 = 1.0 - cflz(i,j,k)
		! 					s2 = 1.0
		! 					ds1 = s2    - s1
		! 				  ds2 = s2**2 - s1**2
		! 			    ds3 = s2**3 - s1**3
		! 			    call ppm(q(i,j,ku-2:ku+2), ql, dq, q6)
		! 			    mfz(i,j,k) = ql * ds1 + 0.5 * dq * ds2 + q6 * (0.5 * ds2 - ds3 / 3.0)
		! 			  else if (cflz(i,j,k) < 0) then
		! 					ku = k 
		! 					s1 = 0
		! 					s2 = -cflz(i,j,k)
		! 				  ds1 = s2    - s1
		! 			    ds2 = s2**2 - s1**2
		! 			    ds3 = s2**3 - s1**3
		! 			    call ppm(q(i,j,ku-2:ku+2), ql, dq, q6)
		! 			    mfz(i,j,k) = -(ql * ds1 + 0.5 * dq * ds2 + q6 * (0.5 * ds2 - ds3 / 3.0))
		! 				else
		! 					mfz(i,j,k) = 0
		! 				end if
		! 			end do
		! 		end do
		!   end do
		end select
		end associate

	end subroutine calc_tracer_mass_flux_vert_ffsl

	subroutine calc_hori_flux_van_leer(tracer, mx, my, mfx, mfy, dt)

		type(tracer_type), intent(in) :: tracer
		real(r8), intent(in) ::   mx(mesh%full_lon_lb:mesh%full_lon_ub, &
														     mesh%full_lat_lb:mesh%full_lat_ub, &
														     mesh%full_lev_lb:mesh%full_lev_ub)
		real(r8), intent(in) ::   my(mesh%full_lon_lb:mesh%full_lon_ub, &
														     mesh%full_lat_lb:mesh%full_lat_ub, &
														     mesh%full_lev_lb:mesh%full_lev_ub)
		real(r8), intent(out) :: mfx(mesh%half_lon_lb:mesh%half_lon_ub, &
														     mesh%full_lat_lb:mesh%full_lat_ub, &
														     mesh%full_lev_lb:mesh%full_lev_ub)
		real(r8), intent(out) :: mfy(mesh%full_lon_lb:mesh%full_lon_ub, &
														     mesh%half_lat_lb:mesh%half_lat_ub, &
														     mesh%full_lev_lb:mesh%full_lev_ub)
		real(r8), intent(in) :: dt
		integer i, j, k, iu, ju, ci
		real(r8) cf, dq

		associate(   u => tracer%u   , v    => tracer%v    ,&
			        cflx => tracer%cflx, cfly => tracer%cfly  )

		! Along x-axis
		do k = mesh%full_lev_ibeg, mesh%full_lev_iend
			do j = mesh%full_lat_ibeg + 1, mesh%full_lat_iend - 1
				do i = mesh%half_lon_ibeg, mesh%half_lon_iend
					ci = int(cflx(i,j,k))
					cf = cflx(i,j,k) - ci
					if (cflx(i,j,k) > 0) then
						iu = i - ci
						dq = slope(mx(iu-1:iu+1,j,k))
						mfx(i,j,k) = (cf * (mx(iu,j,k) + dq * 0.5_r8 * (1 - cf)) + sum(mx(i+1-ci:i,j,k))) / dt * mesh%de_lon(j) 
					else if (cflx(i,j,k) < 0) then
						iu = i - ci + 1
						dq = slope(mx(iu-1:iu+1,j,k))
						mfx(i,j,k) = (cf * (mx(iu,j,k) - dq * 0.5_r8 * (1 + cf)) - sum(mx(i+1:i-ci,j,k))) / dt * mesh%de_lon(j)
					else
						mfx(i,j,k) = 0
					end if
				end do
			end do
		end do
		! Along y-axis
		do k = mesh%full_lev_ibeg, mesh%full_lev_iend
			do j = mesh%half_lat_ibeg, mesh%half_lat_iend
				do i = mesh%full_lon_ibeg, mesh%full_lon_iend
					cf = cfly(i,j,k)
					ju = merge(j, j+1, cf >0)
					dq = slope(my(i,ju-1:ju+1,k))
					mfy(i,j,k) = cf * (my(i,ju,k) + dq * 0.5_r8 * (sign(1.0_r8, cf) - cf)) / dt * mesh%de_lat(j)
				end do
			end do
		end do
		call fill_zonal_halo_lon(mfx, all_halo=.true.)
		end associate

	end subroutine calc_hori_flux_van_leer

	subroutine calc_hori_flux_ppm(tracer, mx, my, mfx, mfy, dt)

		type(tracer_type), intent(in) :: tracer
		real(r8), intent(in) ::   mx(mesh%full_lon_lb:mesh%full_lon_ub, &
														     mesh%full_lat_lb:mesh%full_lat_ub, &
														     mesh%full_lev_lb:mesh%full_lev_ub)
		real(r8), intent(in) ::   my(mesh%full_lon_lb:mesh%full_lon_ub, &
														     mesh%full_lat_lb:mesh%full_lat_ub, &
														     mesh%full_lev_lb:mesh%full_lev_ub)
		real(r8), intent(out) :: mfx(mesh%half_lon_lb:mesh%half_lon_ub, &
														     mesh%full_lat_lb:mesh%full_lat_ub, &
														     mesh%full_lev_lb:mesh%full_lev_ub)
		real(r8), intent(out) :: mfy(mesh%full_lon_lb:mesh%full_lon_ub, &
														     mesh%half_lat_lb:mesh%half_lat_ub, &
														     mesh%full_lev_lb:mesh%full_lev_ub)
		real(r8), intent(in) :: dt
		integer i, j, k, iu, ju, ci
		real(r8) cf, s1, s2, ds1, ds2, ds3, ql, dq, q6

		associate (cflx => tracer%cflx, cfly => tracer%cfly)
		! Along x-axis
			do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				do j = mesh%full_lat_ibeg + 1, mesh%full_lat_iend - 1
					do i = mesh%half_lon_ibeg, mesh%half_lon_iend
						ci = int(cflx(i,j,k))
						cf = cflx(i,j,k) - ci
						if (cflx(i,j,k) > 0) then
							iu = i - ci
							s1 = 1.0 - cf
							s2 = 1.0
						  ds1 = s2    - s1
					    ds2 = s2**2 - s1**2
					    ds3 = s2**3 - s1**3
					    call ppm(mx(iu-2:iu+2,j,k), ql, dq, q6)					
					    mfx(i,j,k) = (sum(mx(i+1-ci:i,j,k)) + ql * ds1 + 0.5 * dq * ds2 + q6 * (0.5 * ds2 - ds3 / 3.0)) / dt * mesh%de_lon(j)
						else if (cflx(i,j,k) < 0) then
							iu = i - ci + 1
							s1 = 0.0
							s2 = -cf
						  ds1 = s2    - s1
					    ds2 = s2**2 - s1**2
					    ds3 = s2**3 - s1**3
					    call ppm(mx(iu-2:iu+2,j,k), ql, dq, q6)
					    mfx(i,j,k) = -(sum(mx(i+1:i-ci,j,k)) + ql * ds1 + 0.5 * dq * ds2 + q6 * (0.5 * ds2 - ds3 / 3.0)) / dt * mesh%de_lon(j)
						else
							mfx(i,j,k) = 0
						end if
					end do
				end do
			end do
			! Along y-axis
			do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				do j = mesh%half_lat_ibeg, mesh%half_lat_iend
					do i = mesh%full_lon_ibeg, mesh%full_lon_iend
						if (cfly(i,j,k) > 0) then
							ju = j
							s1 = 1.0 - cfly(i,j,k)
							s2 = 1.0
							ds1 = s2    - s1
						  ds2 = s2**2 - s1**2
					    ds3 = s2**3 - s1**3
					    call ppm(my(i,ju-2:ju+2,k), ql, dq, q6)
					    mfy(i,j,k) = (ql * ds1 + 0.5 * dq * ds2 + q6 * (0.5 * ds2 - ds3 / 3.0)) / dt * mesh%de_lat(j)
						else if (cfly(i,j,k) < 0) then
							ju = j + 1
							s1 = 0.0
							s2 = -cfly(i,j,k)
							ds1 = s2    - s1
						  ds2 = s2**2 - s1**2
					    ds3 = s2**3 - s1**3
					    call ppm(my(i,ju-2:ju+2,k), ql, dq, q6)
					    mfy(i,j,k) = -(ql * ds1 + 0.5 * dq * ds2 + q6 * (0.5 * ds2 - ds3 / 3.0)) / dt * mesh%de_lat(j)
						else
							mfy(i,j,k) = 0
						end if
					end do
				end do
			end do
			call fill_zonal_halo_lon(mfx, all_halo=.true.)
	  end associate

	end subroutine calc_hori_flux_ppm

	subroutine ppm(f, fl, df, f6)

		real(r8), intent(in) :: f(-2:2)
		real(r8), intent(out) :: fl, df, f6
		real(r8) dfl, dfr, fr

		dfl = slope(f(-2:0))
		df  = slope(f(-1:1))
		dfr = slope(f( 0:2))

		fl = 0.5 * (f(-1) + f(0)) + (dfl - df ) / 6.0
		fr = 0.5 * (f(0)  + f(1)) + (df  - dfr) / 6.0

		fl = f(0) - sign(min(abs(df), abs(fl - f(0))), df)
		fr = f(0) + sign(min(abs(df), abs(fr - f(0))), df)

		f6 = 6 * f(0) - 3 * (fl + fr)
		df = fr - fl

	end subroutine ppm

  pure real(r8) function slope_simple(f) result(res)

  	real(r8), intent(in) :: f(-1:1)
  	real(r8) df

  	df = (f(1) - f(-1)) * 0.5_r8
  	res = df

  end function slope_simple

	pure real(r8) function slope_mono(f) result(res)

	  real(r8), intent(in) :: f(-1:1)
	  real(r8) df, df_min, df_max

	  df = (f(1) - f(-1)) * 0.5_r8
	  df_min = 2 * (f(0) - minval(f))
	  df_max = 2 * (maxval(f) - f(0))
	  res = sign(min(abs(df), df_min, df_max), df) 

	end function slope_mono

	pure real(r8) function slope_pd(f) result(res)

		real(r8), intent(in) :: f(-1:1)
		real(r8) df
		df = (f(1) - f(-1)) * 0.5_r8
		res = sign(min(abs(df), 2 * f(0)), df)

	end function slope_pd

end module ffsl_mod